Nonreciprocal magnetoacoustic waves with out-of-plane phononic angular momenta

Surface acoustic wave (SAW) can carry phononic angular momentum, showing great potential as an energy-efficient way to control magnetism. Still, out-of-plane phononic angular momentum in SAW and its interaction with magnetism remain elusive. Here, we studied the SAW-induced magnetoacoustic waves and spin pumping in Ni-based films on LiNbO3 with selected SAW propagation direction. The crystal inversion asymmetry induces circularly polarized phonons with large out-of-plane angular momenta so that up to 60% of the SAW power attenuates nonreciprocally controlled by the out-of-plane magnetization component. The SAW propagation direction dependence of the nonreciprocity verifies the crystal origin of the phononic angular momentum, and a chiral spin pumping demonstrates that the circular polarization can control the spin current generation efficiency. These results provide an additional degree of freedom for the acoustic control of magnetism and open an avenue for applying circularly polarized phonons.


Theory of SAW-driven FMR
In this section, we derive the theory of SAW-driven FMR from the magnetoelastic coupling between magnetization and strains.In a fully rotational symmetry magnet, which can approximately model the polycrystalline thin films deposited by evaporation, the magnetoelastic free energy reads (18,43) , , , where  =   +   /2 is the linear strain tensor defined on the displacement field u = (ux, uy, uz), b is the magnetoelastic coupling constant, m = (mx, my, mz) is the normalized directional vector of the magnetization.An effective magnetic field is then generated (12,16) where μ0 is the vacuum permeability.For FMR with small precession angle, the effective field can be calculated based on the ground state magnetization mj.
To provide a simple insight on the polarized-phonon-induced FMR, we develop a twodimensional (2D) model with 2D displacement (ux, uy) on a plane (x, y).The displacement and strains are continued at the LiNbO3/Ni interface, and since the thickness of the magnetic film is much thinner than the wavelength, the displacement and strains in the magnetic layer can be approximated by the displacement and strains in the LiNbO3 substrate.The small thickness also means that the interested area is close to the top surface, where the out-of-plane shear strain εxz and εyz should be zero, so that εxz and εyz are small in the Ni layer (14).Meanwhile εyy is zero in a plane wave, and εzz is usually small in LiNbO3 (44).The remaining strains εxx and εxy are only related to the 2D displacement (ux, uy), allowing us to capture the essential physics in the SAWdriven FMR with our 2D model.To account for the states with a tilted magnetization, we define a sphere coordinate system, with magnetization direction m = (sinθcosφ, sinθsinφ, cosθ) (Fig. S1).The h eff components perpendicular to the ground state m are given by (45) As shown in Eq. ( S3) and (S4), both effective field components vanish when θ = 0°, or the magnetization is fully out-of-plane.Therefore, we need to rotate the magnetization from the zaxis.
The linearized Landau-Lifshitz-Gilbert equation takes the form (43,46,47) where (mθ, mφ) = (Δθ, sinθΔφ) is the small derivation of m after Fourier transform in time, γ is the gyromagnetic ratio, α is the Gilbert damping, and F (2) is a self-adjoint two-by-two matrix, relating to the derivation of the magnetic part of the free energy density (modified from Eq. ( 6) in Ref. (45)).For the simplest case of static magnetic field H parallel to the ground state m, 0 1 0 0 1 Solving (Eq.S5) locally and keeping only the linear terms, one can obtain( 16) where ωF = γμ0H.To the leading order, the phonon absorption is given by the work per time per unit volume ΔP done by h eff (16), i.e.
where the star denotes complex conjugation.The second term, either positive or negative, depends on whether the polarization of the effective magnetic field is right-or left-handed circular.Since FMR has a right-hand circular polarization, the right-hand circularly polarized field couples twice as strongly as linear ones, and the left-hand circularly polarized field almost completely decouples.We denote this polarization dependent term ΔPchiral, and replace the effective field with the strains using Eq. ( S3)-(S4) It can be seen that ΔPchiral is maximized when φ = 0° or 180°, which are the chosen experimental conditions in this work.
In a plane wave, the displacement field after Fourier transform in time is given by 0 0 ( , ) ( , ) x y x y z where Lz is given by Eq. (1) in the main text.Combining Eq. (S9)-(S11) and φ = 0°, we obtain Eq. (2) in the main text.

Angular momentum view for the nonreciprocity
We hereby use a semi-classical language to model the nonreciprocity, which more clearly shows the role of the angular momentum.As shown in Ref. (18), considering the predominance strains εxx and εxy only, by introducing m± = mx ± imy, and εx± =εxx ± 2iεxy, the free energy of the magnetoelastic interaction Eq. (S1) can be rewritten as ( ). 2 Note that Ref. (18) defines the in-plane direction perpendicular to SAW propagation as the z axis and discusses an in-plane angular momentum along the same direction, while our z axis is defined along the film normal and we study the out-of-plane angular momentum, so that the above Eq.(S12) has a similar form as Eq.(3) in Ref. (18).The free energy in Eq. ( 12) is also the Hamiltonian of the magnetoelastic interaction, showing that the angular momentum transfer can happen between the magnetic (spin) and the elastic (phonon) systems through the S± and the εx± operators.In the electrical detection of such torque effect, one can analyze the Rayleigh wave-induced FMR in the spin-torque framework (34).Another experiment in YIG/GGG/YIG trilayer also demonstrated the angular momentum transfer between magnons and phonons via the magnetoelastic interaction (48).
We then write down the quantum operator related to phonons and their angular momenta via second quantization.The second quantization of the displacement vector, multiplied by the square root of mass (49), can be written as (SI Eq. (I.5) in Ref. ( where   represents the polarization of the l-th phonon mode, ma is the mass of an atom, with a wave vector k.To study the SAW phonons circularly polarized in the XY plane, we take out the left-and right-handed circular polarized modes with  ± = (1, ±, 0)/√2, wave vector ±k and frequency ω. † l a k and l a k creates and annihilates the general phonons with wave vector k in the l branch, and we focus on the given circularly polarized SAW phonons with creation and annihilation operators † a  and a  , respectively.To describe the surface waves, Eq. ( S13) is based on a two-dimensional plane with N unit cells in total.As we aim at showing the possibility of the angular momentum transfer effect between our SAW phonons and spins, we ignore the contribution from other modes in the following.
By some simple algebra like Eq. (S11), the angular momentum per mass given in Eq. ( 1) in the main text is microscopically equivalent to where ur is the displacement vector of atom at r. Eq. (S14) has a similar form as Eq. ( 1) by Q. Niu et al. (49), while L is an average value for one unit cell, and the displacement is not scaled with the mass of the atom (like in Ref. ( 2)) in our definition.Instead, we only consider the longwavelength phonons (μm scale), so the displacement is the same in one unit cell, allowing us to regard a unit cell as a single "atom" with mass ma.Our expression thus represents the angular momentum per mass.
Replacing u(r) by Eq. (S13), and considering only the XY-circularly polarized SAW with wave vector ±k and frequency ω, one reach (18) . (S15) In other words, Lz is equivalent to the phonon number density difference between the right-and left-circularly polarized phonons.Since the strains are the derivations of displacements, εx± =εxx ± 2iεxy =   ± (  +   ), the operators for εx± are ( 18) As shown in SI Eq. (I.9)-(I.12) in ( 18), these operators raises or lowers the quantum number of the circularly polarized phonon modes, with angular momentum change ħ in each process.To describe the phonon-magnon conversion process, we next rewrite Eq. (S12) with magnon creation and annihilation operators.The equilibrium direction of the magnetization m is tilted by θ from the z axis, requiring a rotation for the coordinate system (Fig. S2).We then rewrite the magnetization in the coordinate system xyz as m = (m x, my, mz) = -(mx' cosθ+ mz'sinθ, my', -mx' sinθ+ mz'cosθ).Here, m' = (mx', my', mz') is the magnetization in the coordinate system x'y'z', which is rotated by θ around the y axis from the coordinate system xyz, and reversed all the axes.The reversion aims at orientating z' along the direction of the spin angular momentum S, which is opposed to m due to the negative charge of electrons.where (S x'(r), Sy'(r), Sz'(r)) is the magnetic spin at r, with spin number S, † b r andb r are the creation and annihilation operators of magnon on site i, respectively.When S is large, there is m' ≈ -(Sx', Sy', Sz')/S.To analyze the magnons with finite wavevectors, the Fourier transformation provides Combining the above expressions, we are now able to take out the contribution of hybridizing magnons and phonons with frequency ω and wavevector ±k in the total magnetoelastic interaction ( ) Note that we simplify (S19) by keeping only the particle number conserving two-particle processes and ignoring all the higher order processes, or the particle number non-conserving processes, as we are working in the linear regime.Also, these particle number conserving twoparticle processes survive after averaging for one period (2π/ω), which cancels out most of the irrelevant processes.  terms convert phonons into magnons, or convert phononic angular momentum (Eq.(S15)) into magnonic angular momentum.Applying the Fermi Golden rule, the phonon to magnon scatter rate can be estimated through the modulus square of the matrix elements, i.e., 2 (2 cos 1) sin As both phonons and magnons carry angular momenta, the scatter rate is also associated with the angular momentum transfer rate.Considering the difference between the scatter rate from the phonons with oppose angular momentum, its θ dependence will be given by (2 cos 1) sin (2 cos 1) sin 4 sin 2 sin , S21) which is the same as derived from the effective field method in S1, displayed in Eq. ( 2) in the main text.
Hence, the SAW nonreciprocity is associated with the difference of the angular momentum transfer rate between the phonons and the magnons, as shown in Fig. S3: With magnetization along -z' and spin along z', when phonons with upward angular momenta are absorbed ( a  ), there is L•S < 0 and it tends to reduce the spin angular momentum, i.e., creating magnons efficiency (Fig. S3A).For phonons with downward angular momenta ( a  ), there is L•S > 0. Although L and S are not fully antiparallel, and it is still possible to create magnons, the scattering efficiency is smaller (Fig. S3B).The efficiency difference is described by the difference in the modulus square of the matrix elements

Sample layout and transmission spectrum
Figure S4A shows the detail of the geometry of the SAW delay line device for the absorption measurement.The separation between port 1 and port 2 IDTs is 300 μm.The IDTs have a width 400 nm and separation 600 nm (wavelength 2 μm).The device used for spin pumping measurement is shown in Fig. S4B.The Pt/Ni heterostructure is patterned into a rectangle with dimension 100 μm × 60 μm. Figure S4C displays the transmission spectrum at zero magnetic field in the ϕk = 150° device for absorption measurement.The magnetic field-dependence data in the main text is collected at the center frequency 1.785 GHz labeled by the black arrow.Center frequency for other devices is taken in a similar way.

Calculation of SAW-driven FMR in Ni film
The magnetization tilting angle θ is calculated using a macrospin model.The magnetic part of the free energy density is 2 0 S (cos cos sin sin cos ) cos , (S22) where H is the applied field, θH is the magnetic field tilting angle in the xz plane, φ is the inplane angle of the magnetization, the saturated magnetization MS ≈ 2.8×10 5 A/m for Ni (51), and anisotropy Kan ≈ μ0MS 2 /2 mainly arising from the demagnetization field.For a given magnetic field H along direction θH, the magnetic ground state can be calculated by minimizing Fm and getting the corresponding θ, with φ always equals to 0°.
The FMR frequency is then calculated from the Landau-Lifshitz equation (52).For the present consideration, one may drop the damping and effective fields from Eq. (S5) to obtain (45) where ωF is the angular frequency of the mode, and the diagonal terms are zero because Taking γ/2π = 28GHz/T, diagonalizing Eq. (S23) with the θ given by minimizing Eq. (S22), we get the FMR eigenfrequency ωF, which is a function of H and θH.With H = 300 mT along θH = 3°, which is the resonance condition in Fig. 2B, Eq. (S23) gives and resonance frequency ωF ≈ 1.8 GHz and θ ≈ 37°.More comprehensive discussion about H, θH and θ, ωF can be found in Ref. (44,51).Now based on the calculated θ and ωF, we can calculate the nonreciprocal power given by Eq. ( 2), which is used for plotting Fig. 2E.The replacement of ωF is based on an approximation that the Ni magnetization under the resonance condition can be modeled by the paramagnetic moments under the same resonance condition and with the same tilting angle θ.We take the excitation frequency ω as the SAW frequency, and the damping constant α = 0.35 for the calculation.The relatively large damping parameter might be caused by the inhomogeneous broadening at low frequency (53,54).Note that the plot in Fig. 2E is in arbitrary unit and does not contain the information of the magnitude of Lz and the magnetoelastic constant b, as we are mainly interested in the evolution of the nonreciprocal absorption with the changing direction and magnitude of the magnetic field.

COMSOL simulation
In COMSOL simulation, we start with the single period simulation with periodic boundary condition, allowing us to analyze the SAW properties in a small unit cell.The unit cell is shown in Fig. S5A, with in-plane dimensions equal to the wavelength λ = 2 μm, and thickness equals to 6λ = 12 μm.A perfect matching layer with thickness 2 μm is set at the bottom to reduce the reflection from the bottom.The mass of the metallic IDTs is ignored.Due to the periodicity, the SAW is standing wave in this simulation.
We first show the admittance as a function of frequency in devices with different orientation in Fig. S5B to find the resonance frequency.The peaks with largest admittance characterize the resonance frequencies 1.827 GHz, 1.781 GHz, 1.835 GHz and 1.937 GHz for ϕk = 0°, 30°, 60° and 90°, respectively.The properties at 120° and 150° are the same as that at 60° and 30°, respectively.
Then we map the strain distributions in these devices.At ϕk = 0° (Fig. S5C), only εxx exists and εxy is zero in the whole device, showing that the transverse and longitudinal modes are decoupled, with only the longitudinal mode couples to the piezoelectricity.At ϕk = 30° or 150° (Fig. S5D), a sizable εxy appears, whose distribution shifts from εxx for approximately a quarter wavelength.Sizable εxy with distribution shifts from εxx also appears along ϕk = 60° or 120° (Fig. S5E), although the shift is smaller than a quarter wavelength.At both 30° and 60°, εxy shifts rightward compared with εxx.When ϕk rotates to 90°, a very small εxy exist, which shifts leftward compared with εxx.Our results are consistent with a recent magnetoacoustics absorption report (55).We then expect large circular polarizations in the propagating waves along ϕk = 30° (150°) and 60° (120°), while 30° (150°) allows a slightly larger circular polarization.Along ϕk = 0°, no circular polarization can exist.A small circular polarization with opposite sign compared to 30° (150°) and 60° (120°) is expected along ϕk = 90°.These correspondence origins from the fact that standing wave is the superposition of propagating waves opposite direction, .It is easy to see that the phase difference between ux0 and uy0, which determines the circular polarization of the propagating waves, is given by the distribution shift of εxx and εxy in the sanding wave.The amplitude of ux0 and uy0, are given by the amplitude of εxx and εxy in the sanding wave.Hence, the angular momentum carried by the propagating waves can be obtained by analyzing the standing waves showing in Fig. S5C-F.We show the emission of the angular-momentum-carrying waves along ϕk = 150° using a larger device structure shown in Fig. 3C.9 IDT fingers are used, with 5 fingers for RF signal and 4 fingers for ground.The wavelength is 2 μm, corresponding to IDT fingers with width 500 nm and separation 500 nm.The z component of the angular momentum in (Eq.1) is calculated as a function of time at the surface of the substrate.The integrated Lz in the AB region and the CD region are used to reflect the angular momentum carried by the + k and -k SAWs, respectively.The simulation along ϕk = 0° is done in a similar way.
We next add the Ni(20)/Ti(15) layers onto the transport part in the above device (Fig S6A) and analyze the strain distribution.In Fig. S6, we show the time-evolution of the strain components in the Ni layer near the LiNbO3 substrate, which is representative for analyzing the magnetoelastic coupling.The time evolutions of the strains in the whole structure are attached in the supplementary animations.As shown in Fig. S6B, the strains εxz, εyz and εyy are negligible.The εyy vanishes as the SAW studied here is a plane wave, and the εxz, εyz are limited by the film thickness (35nm) which is much smaller than the wavelength (2 μm).In addition to the εxx and εxy we have analyzed in the two-dimensional model, a finite εzz appears in the simulation.Taking εzz into account, Eq. (S3) would be modified to 0 0 eff 2 1 sin 2 ( cos sin 2 ), with (S4) remains the same.When φ = 0°, (S24) is equivalent to replacing εxx in Eq. (S3) with a renormalized εxx' = εxx -εzz, and all other following analysis would not be changed.
To evaluate the decay length in the depth direction and the propagation loss, we calculated the acoustic energy flow in the substrate with different SAW propagation direction.As shown in the results (Fig. S6C), while all the SAWs have a decay length comparable to the wavelength in the depth direction, only SAWs with ϕk = 0°(180°) and 90° have good energy confinement at the surface.When ϕk = 30°(150°) or 60°(120°), finite leakage waves propagating into the substrate appears, increasing the propagation loss in these directions.

Piezoelectric effect and phonon circular polarization
We now provide a parsed approximation for a straightforward understanding on how the phononic circular polarization emerges from the piezoelectric effect.The coupling between the displacement components and electric potential in SAWs can be described by the Christoffel equation (56)(57)(58).For a plane SAW at the surface of a piezoelectric material, with propagation direction along x-axis and surface normal along z-axis (Fig. S7A), the Christoffel equation reads where ρ is the mass density of the piezoelectric material, vs is the velocity of the SAW, Γij parameters are related to the elastic, piezoelectric and permittivity matrices, and (Ax, Ay, Az, AU) T describes the wavefunction of the SAW.The displacement and electric potential are given by , (S27) where ω is the circular frequency of the SAW, α is the decay factor along the depth direction (-z direction) of the material.For a surface wave decaying along the depth direction, the imaginary part of α should be positive.The Γij parameters are also related to α as a result of the inhomogeneous distribution of the displacement and the potential.Hence, Eq. ( S25) is an equation about α and vs, and is usually solved in a self-consistent manner, with complex selection rules to distinguish SAW and leaky waves (56,59).Here, to give a qualitative explanation on the phononic angular momentum, we skip this analysis and simply assume α = i, or the studied wave has a delay length equals to the wavelength.
The complex α brings about imaginary parts in Γi4 in Eq. (S28).We then consider a mixing piezoelectric parameter Γmix = 2Im[Γ14 * Γ24], and a normalized Γmix = 2Im[Γ14 * Γ24]/(|Γ14| 2 + |Γ24| 2 ).The relationship between these two parameters and the out-of-plane phononic angular momentum can be understood by considering the following process: when ux, uy are excited by a periodic potential U, their magnitude can be estimated by ux ~ Γ14U, uy ~ Γ24U.Then, the phononic angular momentum ~ Im  3B.The data between 180° and 360° is simply the opposite of the data between 0° and 180°, as displayed in Fig. 1.For most of ϕk, Γmix is nonzero, and large Γmix appears at ϕk between X and Y' axes.At ϕk = 0° (Y' axis), Γmix = 0 because Γ24 = 0, i.e., the transverse deformation uy is decoupled with the electric field and the longitudinal deformation.The normalized Γmix has a maximum value ~ 0.5, consistent with the large circular polarization reflected by the large nonreciprocity in the experiment.The parameter Γmix based on Eq. (S28) also shows a zero value for the Y-cut z-propagating SAW, consistent with a previous study reporting no nonreciprocity (44).Though, the estimation here is simplified by skipping the complicated solving process of the Christoffel equation, aiming at clarifying the existence of the mechanism shown in Fig. 3F, and a quantitative agreement with the experiment shall not be expected based on this analysis.

Fig. S1 .
Fig. S1.Definition of angles characterizing the ground state magnetization orientation.

Fig
Fig.S2.Rotation and inversion of the coordinate system.The coordinate system xyz is transformed to coordinate system x'y'z'.

2 meF
 , as given in Eq. (S21), matching with the effective field helicity mismatch calculation Eq. (S9) and the experimental results.

Fig. S3 .
Fig. S3.The angular momentum view for the magnon-phonon coupling.The phonon-tomagnon scatter rate are different when L•S < 0 (A) and L•S > 0 (B).Note the S fluctuation amplitude ΔS due to magnon generation is exaggerated for a clear view.

Fig. S4 .
Fig. S4.Geometry of the SAW devices.(A) Geometry of the SAW delay line device for the absorption measurement.(B) Geometry of the SAW delay line device for the spin pumping measurement.(C) Transmission spectrum of the SAW delay line device for the absorption measurement with ϕ k = 150°.

Fig. S5 .
Fig. S5.Single-wavelength COMSOL Simulation.(A) Schematic of the unit cell in the single period simulation.(B) Admittance as a function of frequency in devices with different orientation.(C-F) Strain distributions in devices with different orientation.
of εxx and εxy in the sanding wave are given by , [ux * uy] ~ Im[Γ14 * Γ24].To estimate the strength of the circular polarization, we can consider a normalized angular momentum ~ 2Im[ux * uy]/(|ux| 2 + |uy| 2 ) ~ 2Im[Γ14 * Γ24]/(|Γ14| 2 + |Γ24| 2 ).Therefore, Γmix can be used to estimate the strength of the piezoelectric mixing effect that generating the out-of-plane phononic angular momentum, and the normalized Γmix can be used to estimate the strength of the circular polarization.By rotating the piezoelectric matrix [e] of LiNbO3 (60) using the standard Euler and Bond transform matrices (58, 60), we can obtain [e] for the 128° Y-cut LiNbO3 substrate with x-axis along a given direction labeled by ϕk.Fig. S7B plots the corresponding Γmix and normalized Γmix as functions of ϕk on the 128° Y-cut LiNbO3.It can be seen that the Γmix curve can qualitatively reproduce the tendency of the experimental ΔPn curve in Fig.

Fig. S7 .
Fig. S7.Piezoelectric strain-mixing effect.(A) Schematic of the coordinate system for the Christoffel equation of the SAW.(B) Propagation direction dependence of the mixing piezoelectric parameters.